Biological analysis
2019/04/08
1 Introduction
We already demultiplexed and assigned cells back to their original condition, filtered low-quality cells, normalized UMI counts and classified cells by cell types (CD4 T , Cytotoxic (CD8 T + NK), monocytes and B cells). Thus, we dispose of a processed dataset to start answering the following questions:
- Does the time required to freeze a sample bias a cell’s transcriptome (i.e. introduces technical artifacts)?
- If (1) is true, can we correct it by adjusting the experimental conditions (i.e. changing the temperature)?
1.1 Package loading
library(Matrix)
library(stringr)
library(psych)
library(kmed)
library(pheatmap)
library(fitdistrplus)
library(ggpubr)
library(SingleCellExperiment)
library(scater)
library(scran)
library(SC3)
library(Seurat)
library(purrr)
library(grid)
library(gridExtra)
library(gridGraphics)
library(Seurat)
library(cowplot)
library(tidyverse)1.2 Source function definitions
source(file = "bin/utils.R")1.3 Load SingleCellExperiment
pbmc <- readRDS("results/R_objects/10X_SingleCellExperiment_clustered.RDS")2 Effect of time on gene expression
First, we aim to elucidate whether the time that a blood sample spends at room temperature before cryopreservation introduces technical noise to the gene expression matrix. We dispose of data from two donors (male and female), which we will use as biological replicates. For each of them, the same blood sample was kept at room temperature for several times: 0h, 2h, 8h, 24h, 48h. Our hypothesis is that, if this variable introduces artifacts, cells will cluster by time. Thus, our first approach will be to perform non-linear dimensionality reduction (tSNE) to visualize such an effect.
date <- Sys.Date()
# Separate SCE into the two replicates: male and female
pbmc_l <- list(male = pbmc[, pbmc$sex == "male"], female = pbmc[, pbmc$sex == "female"])
# Filter out 4ºC cells, as we will focus on them later
pbmc_l <- map(pbmc_l, ~ .[, .$temperature != "4ºC"])
# Find hypervariable genes
pbmc_l <- map(pbmc_l, function(x) {
fit_var <- trendVar(x, use.spikes = FALSE)
decomp_var <- decomposeVar(x, fit_var)
top_hvgs <- order(decomp_var$bio, decreasing = TRUE)
top_20_pct_hvgs <- top_hvgs[1:(0.2 * length(top_hvgs))]
x[top_20_pct_hvgs, ]
})
# Plot tSNE
time_points <- c("0h", "2h", "8h", "24h", "48h")
pbmc_l$female$time <- factor(pbmc_l$female$time, levels = time_points)
pbmc_l$male$time <- factor(pbmc_l$male$time, levels = time_points)
palette <- c("#999999", "#92e8df", "#632c63", "#e4624e", "#c0e212")
set.seed(1)
tsnes_sex <- map(c("male", "female"), function(sex) {
plot_tsne(
pbmc_l[[sex]],
exprs_values = "logcounts",
color_by = "time",
colors = palette,
point_size = 0.5,
point_alpha = 1,
title = str_to_title(sex)
)
})
tsnes_sex <- map(tsnes_sex, ~ . +
guides(colour = guide_legend(override.aes = list(size = 4))))
tsnes1 <- ggarrange(
plotlist = tsnes_sex,
nrow = 2,
ncol = 1,
common.legend = TRUE
)
# Save plot
ggsave(
plot = tsnes1,
filename = str_c("results/plots/", date, "_tSNEs_male&female.pdf"),
device = "pdf",
width = 12,
height = 6.5
)
saveRDS(object = tsnes_sex[[2]], "results/R_objects/tsne_time_points_female_gg.rds")
tsnes_sex## [[1]]
##
## [[2]]
Indeed, we see how time until cryopreservation has an effect on gene-expression. Specifically, we see a clear gradient from 0h-2h to 48h.
Let us assess whether time is explained by any of the principal components (PCs) for each cell type:
# Separate by cell types
pbmc_fil <- pbmc[, pbmc$temperature != "4ºC"]
pbmc_types <- list(
"CD4 T" = pbmc_fil[, pbmc_fil$cell_type == "CD4 T"],
"CD8 T" = pbmc_fil[, pbmc_fil$cell_type == "CD8 T"],
"NK" = pbmc_fil[, pbmc_fil$cell_type == "NK"],
"Monocyte" = pbmc_fil[, pbmc_fil$cell_type == "Monocyte"],
"B" = pbmc_fil[, pbmc_fil$cell_type == "B"]
)
# Run PCA
pbmc_types <- pbmc_types %>%
map(find_var_genes) %>%
map(~ runPCA(., ntop = nrow(.)))
# Plot PC1 vs time
pc1_time_df <- map(pbmc_types, function(sce) {
reducedDim(sce, "PCA") %>%
as.data.frame() %>%
set_names(c("PC1", "PC2")) %>%
dplyr::mutate(time = sce$time)
})
pc1_time_df <- bind_rows(pc1_time_df, .id = "cell_type")
palette2 <- c("#c20a35", "#aa2edc", "green2", "#bbaa2a", "#71bdd0")
pc1_time_gg <- pc1_time_df %>%
dplyr::mutate(time = factor(time, c("0h", "2h", "8h", "24h", "48h"))) %>%
dplyr::mutate(cell_type = factor(cell_type, c("CD4 T", "CD8 T", "NK", "Monocyte", "B"))) %>%
filter(PC1 > -10, PC1 < 10) %>%
ggplot(aes(time, PC1, fill = cell_type)) +
geom_boxplot(outlier.shape = NA) +
scale_fill_manual(values = palette2) +
labs(x = "", fill = "") +
theme_classic()
saveRDS(object = pc1_time_gg, file = "results/R_objects/pc1_time_gg.rds")
pc1_time_gg
ggsave(
plot = pc1_time_gg,
filename = str_c("results/plots/", date, "_pc1_vs_time.pdf"),
width = 9,
height = 7
)
# #version 2#####
# pc1_time_gg2 <- pc1_time_df %>%
# dplyr::mutate(time = factor(time, c("0h", "2h", "8h", "24h", "48h"))) %>%
# dplyr::mutate(cell_type = factor(cell_type, c("CD4 T", "CD8 T", "NK", "Monocyte", "B"))) %>%
# filter(PC1 > -10, PC1 < 10) %>%
# ggplot(aes(time, PC1)) +
# geom_boxplot(outlier.shape = NA, fill = "grey") +
# facet_grid(. ~ cell_type) +
# labs(x = "") +
# theme_bw()
# pc1_time_gg2
# ggsave(
# plot = pc1_time_gg2,
# filename = str_c("results/plots/", date, "_pc1_vs_time2.pdf"),
# width = 11,
# height = 5
# )Interestingly, for most PBMC types, processing time represents the major source of variability (as it explains the variability in PC1). This is not as obvious for B cells, likely because there are too few.
2.1 Zoom-in into each cell type
We are also interested in whether the aforementioned effect is cell-type specific or whether is ubiquitously found in all cell types:
# For each sex and cell_type, run tSNE to zoom-in into specific populations
set.seed(2)
cell_types <- levels(pbmc$cell_type)
zoom_l <- map(c("male", "female"), function(x) {
tsne_df <- data.frame()
for (cell_type in cell_types) {
curr_pbmc <- pbmc_l[[x]][, pbmc_l[[x]]$cell_type == cell_type]
curr_pbmc <- find_var_genes(curr_pbmc)
curr_pbmc <- runTSNE(object = curr_pbmc, exprs_values = "logcounts")
curr_df <- reducedDim(curr_pbmc, "TSNE") %>%
as.data.frame() %>%
set_names(c("TSNE1", "TSNE2")) %>%
mutate(time = colData(curr_pbmc)[["time"]],
cell_type = rep(cell_type, ncol(curr_pbmc)))
tsne_df <- rbind(tsne_df, curr_df)
}
tsne_df
})
# Join male and female into same df
# Relevel categorical variables for visualization
zoom_df <- zoom_l %>%
set_names(c("male", "female")) %>%
bind_rows(.id = "sex") %>%
mutate(sex = factor(sex, levels = c("male", "female")),
cell_type = factor(cell_type, levels = cell_types),
time = factor(time, levels = c("0h", "2h", "8h", "24h", "48h")))
# Plot tSNE
zoom_tsne <- zoom_df %>%
ggplot(aes(TSNE1, TSNE2, color = time)) +
geom_point(size = 0.2, alpha = 1) +
scale_color_manual(values = palette) +
labs(color = "") +
facet_grid(sex ~ cell_type) +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
panel.grid = element_blank(),
strip.text = element_text(size = 12),
legend.text = element_text(size = 12)) +
guides(colour = guide_legend(override.aes = list(size = 4)))
# Save tSNE
ggsave(
plot = zoom_tsne,
filename = str_c("results/plots/", date, "_tSNEs_zoom.pdf"),
device = "pdf",
height = 7.5,
width = 15
)
saveRDS(object = zoom_tsne, "results/R_objects/tsne_time_points_by_cell_types_gg.rds")
zoom_tsne
We can clearly visualize that not only is the effect consistent across the major cell types, but also in subtypes that contain fewer cells.
3 Room temperature vs 4ºC
Finally, we want to inspect if maintaining cells at 4ºC prevents the bias from appearing:
# Separate SCE into the two replicates: male and female
pbmc_temp <- list(
male = pbmc[, pbmc$sex == "male"],
female = pbmc[, pbmc$sex == "female"]
)
# Filter out 2h
pbmc_temp <- map(pbmc_temp, ~ .[, .$time != "2h"])
# Find hypervariable genes
pbmc_temp <- map(pbmc_temp, function(x) {
fit_var <- trendVar(x, use.spikes = FALSE)
decomp_var <- decomposeVar(x, fit_var)
top_hvgs <- order(decomp_var$bio, decreasing = TRUE)
top_20_pct_hvgs <- top_hvgs[1:(0.2 * length(top_hvgs))]
x[top_20_pct_hvgs, ]
})
# Plot tSNE
temps <- c("gold", "room temperature", "4ºC")
pbmc_temp$female$temperature <- factor(pbmc_temp$female$temperature, temps)
pbmc_temp$male$temperature <- factor(pbmc_temp$male$temperature, temps)
tsnes_temp <- map(c("male", "female"), function(sex) {
plot_tsne(
pbmc_temp[[sex]],
exprs_values = "logcounts",
color_by = "temperature",
colors = c("#999999", "#ed9121", "#a5cded"),
point_size = 0.5,
point_alpha = 0.9,
title = str_to_title(sex)
)
})
tsnes_temp <- map(tsnes_temp, ~ . +
guides(colour = guide_legend(override.aes = list(size = 4))))
tsnes_temp2 <- ggarrange(
plotlist = tsnes_temp,
nrow = 1,
ncol = 2,
common.legend = TRUE
)
# Save plot
ggsave(
plot = tsnes_temp2,
filename = str_c("results/plots/", date, "_tSNEs_tempearature.pdf"),
device = "pdf",
height = 7,
width = 12
)
saveRDS(tsnes_temp[[2]], file = "results/R_objects/tsne_temperature_female.rds")
tsnes_temp## [[1]]
##
## [[2]]
As we can see, cells that were kept at 4ºC before cryopreservation cluster with the “gold-standard” cells (i.e. those that were cryopreserved right after extraction). Hence, we can conclude that cells should be kept at the fridge before processing. Note however, that we have very few 4ºC cells for the male donor, so we cannot be absolutely certain that this is consistent across individuals. Moreover, although the ‘mixedness’ is much better, 4ºC cells still seem to mix preferably between them than with 0h cells.
Let us assess whether there are any differentially expressed genes between 4ºC and 0h:
# Filter pbmc to retain only 4ºC and 0h cells from the female donor
pbmc_temp_f <- pbmc[, pbmc$sex == "female" & pbmc$temperature != "room temperature"]
# Split by cell type
pbmc_temp_l <-list(
"CD4 T" = pbmc_temp_f[, pbmc_temp_f$cell_type == "CD4 T"],
"CD8 T" = pbmc_temp_f[, pbmc_temp_f$cell_type == "CD8 T"],
"NK" = pbmc_temp_f[, pbmc_temp_f$cell_type == "NK"],
"Monocyte" = pbmc_temp_f[, pbmc_temp_f$cell_type == "Monocyte"],
"B" = pbmc_temp_f[, pbmc_temp_f$cell_type == "B"]
)
# Convert to Seurat objects
pbmc_temp_l <- map(pbmc_temp_l, Convert, to = "seurat")
pbmc_temp_l <- map(pbmc_temp_l, ScaleData)## NormalizeData has not been run, therefore ScaleData is running on non-normalized values. Recommended workflow is to run NormalizeData first.
## ScaleData is running on non-normalized values. Recommended workflow is to run NormalizeData first.
## NormalizeData has not been run, therefore ScaleData is running on non-normalized values. Recommended workflow is to run NormalizeData first.
## ScaleData is running on non-normalized values. Recommended workflow is to run NormalizeData first.
## NormalizeData has not been run, therefore ScaleData is running on non-normalized values. Recommended workflow is to run NormalizeData first.
## ScaleData is running on non-normalized values. Recommended workflow is to run NormalizeData first.
## NormalizeData has not been run, therefore ScaleData is running on non-normalized values. Recommended workflow is to run NormalizeData first.
## ScaleData is running on non-normalized values. Recommended workflow is to run NormalizeData first.
## NormalizeData has not been run, therefore ScaleData is running on non-normalized values. Recommended workflow is to run NormalizeData first.
## ScaleData is running on non-normalized values. Recommended workflow is to run NormalizeData first.pbmc_temp_l <- map(pbmc_temp_l, SetAllIdent, id = "temperature")
output_dea <- map(pbmc_temp_l, function(seurat) {
dea_output <- FindMarkers(
seurat,
ident.1 = "4ºC",
ident.2 = "gold",
test.use = "MAST"
)
})
# Heatmap DEG
output_dea_sig <- map(output_dea, function(df) {
df <- df %>%
rownames_to_column(var = "gene") %>%
filter(p_val_adj < 0.01 & avg_logFC > 1)
rownames(df) <- df$gene
df
})
all_sig_genes <- output_dea_sig %>%
map("gene") %>%
unlist() %>%
unique()
p_adj_list <- map(all_sig_genes, function(gene) {
p_adj_vec <- map_dbl(cell_types, ~ output_dea_sig[[.]][gene, "p_val_adj"] )
p_adj_vec[is.na(p_adj_vec)] <- 1
p_adj_vec
})
names(p_adj_list) <- all_sig_genes
p_adj_mat <- p_adj_list %>%
bind_rows() %>%
t() %>%
as.matrix()
colnames(p_adj_mat) <- cell_types
p_adj_mat <- -1 * log10(p_adj_mat)
cols <- colorRampPalette(c("gray99", "brown2"))(30)
p_adj_mat <- p_adj_mat[, c(1:3, 5, 4)]
heat_4vs0 <- pheatmap(p_adj_mat, cluster_rows = FALSE, angle_col = 45, cluster_cols = FALSE, color = cols)
# Save
pdf(
file = str_c("results/plots/", date, "_heatmap_4ºC_vs_0h.pdf"),
height = 6,
width = 3
)
print(heat_4vs0)
dev.off()## quartz_off_screen
## 2pdf(
file = str_c("doc/figures/R/", date, "_heatmap_4ºC_vs_0h.pdf"),
height = 6.9,
width = 3.3
)
print(heat_4vs0)
dev.off()## quartz_off_screen
## 2Although very few genes were differentially expressed, they are highly biologically meaningful. For instance, JUN, a gene essential for stress adaptation (as explained in this paper), was differentially expressed in all PBMC types.
4 Final figure
# Figure effect time at RT and 4ºC
tsne_time <- tsnes_sex[[2]] +
ggtitle(NULL) +
theme(axis.line = element_blank(),
axis.ticks = element_blank(),
axis.text = element_blank(),
axis.title = element_blank(),
panel.border = element_blank(),
legend.position = "top")
tsne_temp <- tsnes_temp[[2]] +
ggtitle(NULL) +
theme(axis.line = element_blank(),
axis.ticks = element_blank(),
axis.text = element_blank(),
axis.title = element_blank(),
panel.border = element_blank(),
legend.position = "top") +
guides(colour = guide_legend(override.aes = list(size = 4)))
tsne_zoom <- zoom_df %>%
filter(sex == "female") %>%
ggplot(aes(TSNE1, TSNE2, color = time)) +
geom_point(size = 2, alpha = 1) +
scale_color_manual(values = palette) +
labs(color = "") +
facet_grid(. ~ cell_type) +
theme_bw() +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
panel.grid = element_blank(),
strip.text = element_text(size = 12),
legend.position = "none") +
guides(colour = guide_legend(override.aes = list(size = 4)))
tsne_time_temp <- ggarrange(plotlist = list(tsne_time, tsne_temp), ncol = 2, nrow = 1)
effect_time_figure <- plot_grid(
tsne_time_temp,
NULL,
tsne_zoom,
nrow = 3,
ncol = 1,
rel_heights = c(1.2, 0.1, 0.8)
)
ggsave(
filename = str_c("results/plots/", date, "_effect_time_figure.pdf"),
plot = effect_time_figure,
width = 12,
height = 12
)
effect_time_figure
4.1 Supplementary figure
leg <- as_ggplot(get_legend(tsnes_sex[[1]]))
tsne_male <- tsnes_sex[[1]] +
theme(plot.title = element_blank(),
legend.position = "none",
plot.background = element_blank(),
panel.border = element_blank(),
axis.title = element_blank(),
axis.ticks = element_blank(),
axis.text = element_blank())
zoom_tsne <- zoom_tsne +
theme_classic() +
theme(legend.position = "none",
axis.text = element_text(size = 7),
axis.title = element_text(size = 10)) +
labs(x = "tSNE1", y = "tSNE2")
reproducible_donors_gg <- plot_grid(
tsne_male,
NULL,
zoom_tsne,
nrow = 1,
ncol = 3,
rel_widths = c(0.35, 0.02, 0.75),
labels = c("a", "", "b")
)
ggsave(
filename = str_c("doc/figures/legends/", date, "_legend_tsne_time.pdf"),
plot = leg,
width = 1,
height = 3
)
ggsave(
filename = str_c("doc/figures/R/", date, "_reproducible_donors.pdf"),
plot = reproducible_donors_gg,
width = 18.5,
height = 7,
units = "cm"
)5 Session Info
sessionInfo()## R version 3.5.1 (2018-07-02)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] grid parallel stats4 stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] forcats_0.4.0 dplyr_0.8.0.1 readr_1.3.1 tidyr_0.8.3 tibble_2.1.1 tidyverse_1.2.1 gridGraphics_0.3-0 gridExtra_2.3 purrr_0.3.2 Seurat_2.3.4 cowplot_0.9.4 SC3_1.10.1 scran_1.10.2 scater_1.10.1 SingleCellExperiment_1.4.1 SummarizedExperiment_1.12.0 DelayedArray_0.8.0 BiocParallel_1.16.6 matrixStats_0.54.0 Biobase_2.42.0 GenomicRanges_1.34.0 GenomeInfoDb_1.18.2 IRanges_2.16.0 S4Vectors_0.20.1 BiocGenerics_0.28.0 ggpubr_0.2 magrittr_1.5 ggplot2_3.1.0 fitdistrplus_1.0-14 npsurv_0.4-0 lsei_1.2-0 survival_2.44-1.1 MASS_7.3-51.3 pheatmap_1.0.12 kmed_0.2.0 psych_1.8.12 stringr_1.4.0 Matrix_1.2-17 BiocStyle_2.10.0
##
## loaded via a namespace (and not attached):
## [1] reticulate_1.11.1 R.utils_2.8.0 tidyselect_0.2.5 htmlwidgets_1.3 trimcluster_0.1-2.1 Rtsne_0.15 munsell_0.5.0 codetools_0.2-16 ica_1.0-2 statmod_1.4.30 withr_2.1.2 colorspace_1.4-1 knitr_1.22 rstudioapi_0.10 ROCR_1.0-7 robustbase_0.93-4 dtw_1.20-1 gbRd_0.4-11 labeling_0.3 Rdpack_0.10-1 lars_1.2 GenomeInfoDbData_1.2.0 mnormt_1.5-5 bit64_0.9-7 rhdf5_2.26.2 generics_0.0.2 xfun_0.5 diptest_0.75-7 R6_2.4.0 doParallel_1.0.14 ggbeeswarm_0.6.0 locfit_1.5-9.1 hdf5r_1.0.1 flexmix_2.3-15 bitops_1.0-6 assertthat_0.2.1 promises_1.0.1 SDMTools_1.1-221 scales_1.0.0 nnet_7.3-12 beeswarm_0.2.3 gtable_0.3.0 rlang_0.3.3 splines_3.5.1 lazyeval_0.2.2 acepack_1.4.1 broom_0.5.1
## [48] checkmate_1.9.1 abind_1.4-5 modelr_0.1.4 BiocManager_1.30.4 yaml_2.2.0 reshape2_1.4.3 backports_1.1.3 httpuv_1.5.0 Hmisc_4.2-0 tools_3.5.1 bookdown_0.9 gplots_3.0.1.1 RColorBrewer_1.1-2 proxy_0.4-23 dynamicTreeCut_1.63-1 ggridges_0.5.1 Rcpp_1.0.1 plyr_1.8.4 progress_1.2.0 base64enc_0.1-3 zlibbioc_1.28.0 RCurl_1.95-4.12 prettyunits_1.0.2 rpart_4.1-13 pbapply_1.4-0 viridis_0.5.1 zoo_1.8-5 haven_2.1.0 cluster_2.0.7-1 data.table_1.12.0 lmtest_0.9-36 RANN_2.6.1 mvtnorm_1.0-10 hms_0.4.2 mime_0.6 evaluate_0.13 xtable_1.8-3 readxl_1.3.1 mclust_5.4.3 compiler_3.5.1 KernSmooth_2.23-15 crayon_1.3.4 R.oo_1.22.0 htmltools_0.3.6 segmented_0.5-3.0 pcaPP_1.9-73 later_0.8.0
## [95] Formula_1.2-3 snow_0.4-3 rrcov_1.4-7 lubridate_1.7.4 WriteXLS_4.1.0 fpc_2.1-11.1 MAST_1.8.2 cli_1.1.0 R.methodsS3_1.7.1 gdata_2.18.0 metap_1.1 igraph_1.2.4 pkgconfig_2.0.2 registry_0.5-1 foreign_0.8-71 xml2_1.2.0 foreach_1.4.4 vipor_0.4.5 rngtools_1.3.1 pkgmaker_0.27 XVector_0.22.0 rvest_0.3.2 bibtex_0.4.2 doRNG_1.7.1 digest_0.6.18 tsne_0.1-3 cellranger_1.1.0 rmarkdown_1.12 htmlTable_1.13.1 edgeR_3.24.3 DelayedMatrixStats_1.4.0 kernlab_0.9-27 shiny_1.2.0 gtools_3.8.1 modeltools_0.2-22 jsonlite_1.6 nlme_3.1-137 Rhdf5lib_1.4.3 BiocNeighbors_1.0.0 viridisLite_0.3.0 limma_3.38.3 pillar_1.3.1 lattice_0.20-38 httr_1.4.0 DEoptimR_1.0-8 glue_1.3.1 png_0.1-7
## [142] prabclus_2.2-7 iterators_1.0.10 bit_1.1-14 mixtools_1.1.0 class_7.3-15 stringi_1.4.3 HDF5Array_1.10.1 doSNOW_1.0.16 latticeExtra_0.6-28 caTools_1.17.1.2 irlba_2.3.3 e1071_1.7-1 ape_5.3